{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Graphical Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import fetch_mldata\n",
    "from prml import bayesnet as bn\n",
    "\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = bn.discrete([0.1, 0.9])\n",
    "f = bn.discrete([0.1, 0.9])\n",
    "\n",
    "g = bn.discrete([[[0.9, 0.8], [0.8, 0.2]], [[0.1, 0.2], [0.2, 0.8]]], b, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.1 0.9])\n",
      "f: DiscreteVariable(proba=[0.1 0.9])\n",
      "g: DiscreteVariable(proba=[0.315 0.685])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "f: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(observed=[1. 0.])\n",
      "f: DiscreteVariable(proba=[0.11111111 0.88888889])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.3 Illustration: Image de-noising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ada5ce4c18>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAC0tJREFUeJzt3UGoZGeZxvH/M1E3MYsOIU0Tk4kjYTYu4tC4UYaehZJx03GRwaxaZtEuJqA7g5sERJBBndkJGWzsgTESiJomDBODOBNXIZ0gpmNPTJCe2KbpJvTCZCWadxb3tFw7996qW1WnTt1+/z8oqurcunXerttPfd853znnS1UhqZ+/mLoASdMw/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmnrfOleWxMMJpZFVVeZ53VItf5L7krya5PUkDy/zXpLWK4se25/kJuBXwKeAi8ALwINV9cs9fseWXxrZOlr+jwOvV9Wvq+r3wPeB40u8n6Q1Wib8dwC/2fb84rDszyQ5meRskrNLrEvSii2zw2+nrsV7uvVV9RjwGNjtlzbJMi3/ReDObc8/BLy5XDmS1mWZ8L8A3JPkw0k+AHwOOLOasiSNbeFuf1X9IclDwDPATcCpqnplZZVJGtXCQ30Lrcxtfml0aznIR9LBZfilpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU2udolv9jHl16GSui9RqF7b8UlOGX2rK8EtNGX6pKcMvNWX4paYMv9TUUuP8SS4AbwN/BP5QVUdXUZQOjnXO8rzfdXscwN5WcZDP31XVWyt4H0lrZLdfamrZ8Bfw4yQvJjm5ioIkrcey3f5PVNWbSW4Hnk3yv1X13PYXDF8KfjFIGyar2mGT5FHgnar6xh6vmW7vkEYx5Q6/Wbru8Kuquf7hC3f7k9yc5JZrj4FPA+cWfT9J67VMt/8w8MPh2/V9wPeq6r9WUpWk0a2s2z/Xyuz2Hzib3K2fxW7/3hzqk5oy/FJThl9qyvBLTRl+qSnDLzXlpbub2+ShvFlDdbNq3+vnXYcBt7Pll5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmHOe/AWzyWP2m8rLftvxSW4ZfasrwS00Zfqkpwy81Zfilpgy/1JTj/FrKsufcH9R13whs+aWmDL/UlOGXmjL8UlOGX2rK8EtNGX6pqZnhT3IqyZUk57YtuzXJs0leG+4PjVtmb1W1521MSfa8Lfv7y7y3ljNPy/9d4L7rlj0M/KSq7gF+MjyXdIDMDH9VPQdcvW7xceD08Pg0cP+K65I0skW3+Q9X1SWA4f721ZUkaR1GP7Y/yUng5NjrkbQ/i7b8l5McARjur+z2wqp6rKqOVtXRBdclaQSLhv8McGJ4fAJ4ajXlSFqXzHEJ48eBY8BtwGXgEeBHwBPAXcAbwANVdf1OwZ3ey3MsFzDlqakHechtmc/tgP+75yp+ZvhXyfDvzHCPY8prCUxp3vB7hJ/UlOGXmjL8UlOGX2rK8EtNGX6pKS/dvQbrOO1W2i9bfqkpwy81Zfilpgy/1JThl5oy/FJThl9qynH+A8BxfI3Bll9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmnKcfwWmvPS2tChbfqkpwy81Zfilpgy/1JThl5oy/FJThl9qamb4k5xKciXJuW3LHk3y2yQ/H26fGbfMG1uSPW/SGOZp+b8L3LfD8n+pqnuH23+utixJY5sZ/qp6Dri6hlokrdEy2/wPJfnFsFlwaGUVSVqLRcP/beAjwL3AJeCbu70wyckkZ5OcXXBdkkaQeU5KSXI38HRVfXQ/P9vhtTfkGTDLntjjTr1xjHnC1Sb/zapqruIWavmTHNn29LPAud1eK2kzzTylN8njwDHgtiQXgUeAY0nuBQq4AHxhxBoljWCubv/KVma3f0eb3IXcZF279bOM2u2XdPAZfqkpwy81Zfilpgy/1JThl5ry0t3aWA7ljcuWX2rK8EtNGX6pKcMvNWX4paYMv9SU4Zeacpxfkxn7dHLH8vdmyy81Zfilpgy/1JThl5oy/FJThl9qyvBLTTnOr1F5Tv7msuWXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paZmhj/JnUl+muR8kleSfHFYfmuSZ5O8NtwfGr/cG1NV7Xmbct3L3paRZM+blpNZf6AkR4AjVfVSkluAF4H7gc8DV6vq60keBg5V1ZdnvNe4/5MncpAvSjF27csw4Iupqrk+uJktf1VdqqqXhsdvA+eBO4DjwOnhZafZ+kKQdEDsa5s/yd3Ax4DngcNVdQm2viCA21ddnKTxzH1sf5IPAk8CX6qq383bJUtyEji5WHmSxjJzmx8gyfuBp4Fnqupbw7JXgWNVdWnYL/DfVfXXM95nczcwl+A2/zjc5l/Myrb5s/UX+A5w/lrwB2eAE8PjE8BT+y1S0nTm2dv/SeBnwMvAu8Pir7C13f8EcBfwBvBAVV2d8V6b28wsYZNbz01myz6OeVv+ubr9q2L4tZ3hH8fKuv2SbkyGX2rK8EtNGX6pKcMvNWX4paa8dLeW4nDdwWXLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNOc6/ArPGujf5lF/H6fuy5ZeaMvxSU4ZfasrwS00Zfqkpwy81ZfilphznXwPH0rWJbPmlpgy/1JThl5oy/FJThl9qyvBLTRl+qamZ4U9yZ5KfJjmf5JUkXxyWP5rkt0l+Ptw+M365klYlsy40keQIcKSqXkpyC/AicD/wD8A7VfWNuVeWbO5VLaQbRFXNdVTZzCP8quoScGl4/HaS88Ady5UnaWr72uZPcjfwMeD5YdFDSX6R5FSSQ7v8zskkZ5OcXapSSSs1s9v/pxcmHwT+B/haVf0gyWHgLaCAr7K1afCPM97Dbr80snm7/XOFP8n7gaeBZ6rqWzv8/G7g6ar66Iz3MfzSyOYN/zx7+wN8Bzi/PfjDjsBrPguc22+RkqYzz97+TwI/A14G3h0WfwV4ELiXrW7/BeALw87Bvd7Lll8a2Uq7/ati+KXxrazbL+nGZPilpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2pq3VN0vwX837bntw3LNtGm1rapdYG1LWqVtf3lvC9c6/n871l5craqjk5WwB42tbZNrQusbVFT1Wa3X2rK8EtNTR3+xyZe/142tbZNrQusbVGT1DbpNr+k6Uzd8kuayCThT3JfkleTvJ7k4Slq2E2SC0leHmYennSKsWEatCtJzm1bdmuSZ5O8NtzvOE3aRLVtxMzNe8wsPelnt2kzXq+925/kJuBXwKeAi8ALwINV9cu1FrKLJBeAo1U1+Zhwkr8F3gH+/dpsSEn+GbhaVV8fvjgPVdWXN6S2R9nnzM0j1bbbzNKfZ8LPbpUzXq/CFC3/x4HXq+rXVfV74PvA8Qnq2HhV9Rxw9brFx4HTw+PTbP3nWbtdatsIVXWpql4aHr8NXJtZetLPbo+6JjFF+O8AfrPt+UU2a8rvAn6c5MUkJ6cuZgeHr82MNNzfPnE915s5c/M6XTez9MZ8dovMeL1qU4R/p9lENmnI4RNV9TfA3wP/NHRvNZ9vAx9haxq3S8A3pyxmmFn6SeBLVfW7KWvZboe6Jvncpgj/ReDObc8/BLw5QR07qqo3h/srwA/Z2kzZJJevTZI63F+ZuJ4/qarLVfXHqnoX+Dcm/OyGmaWfBP6jqn4wLJ78s9uprqk+tynC/wJwT5IPJ/kA8DngzAR1vEeSm4cdMSS5Gfg0mzf78BngxPD4BPDUhLX8mU2ZuXm3maWZ+LPbtBmvJznIZxjK+FfgJuBUVX1t7UXsIMlfsdXaw9YZj9+bsrYkjwPH2Drr6zLwCPAj4AngLuAN4IGqWvuOt11qO8Y+Z24eqbbdZpZ+ngk/u1XOeL2SejzCT+rJI/ykpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzX1//kg+AXbkzT0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mnist = fetch_mldata(\"MNIST original\")\n",
    "x = mnist.data[0]\n",
    "binarized_img = (x > 127).astype(np.int).reshape(28, 28)\n",
    "plt.imshow(binarized_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ada5d84898>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADMBJREFUeJzt3U+IpHedx/H3d6NeYg4JMbNDTHZcCV5yiNuNF2UZD0oUYeIhwZxGhG0PBtabIZcEFiEs6q4nYXYdHEGjgagZwrJRRI2nkJ6wmOj4J8hsnE0zo4xgchLNdw/9jLST7qrqep6nnuep7/sFTVfX1NTzrafqU8+vnu/z1C8yE0n1/M3QBUgahuGXijL8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1TUm1a5sIhYy8MJNzY2Zv77uXPnBlt+38ueqqGfsz5lZixyu2hzeG9E3A18EbgO+M/MfHTO7dcy/PPWYcRCz0Uvy+972VM19HPWp97DHxHXAb8EPgBcBJ4D7s/Mn834P4Z/xcuf8ou4T0M/Z31aNPxtPvO/B3gpM3+dmX8EvgGcaHF/klaoTfhvBX6z5++LzXV/JSK2ImI7IrZbLEtSx9rs8NtvaPGGsVRmngJOwfoO+6UparPlvwjctufvtwOvtCtH0qq0Cf9zwB0R8Y6IeAvwMeBsN2VJ6tvSw/7M/FNEPAA8zW6r73Rm/rSzyiZk6D3Ds5a/znu126j6uPdq1ec/9ML8zL9yhr+eVbT6JE2Y4ZeKMvxSUYZfKsrwS0UZfqmolZ7PP6S2Lc0+W2J9tuNs5a2fWa+Xzc3Nhe/HLb9UlOGXijL8UlGGXyrK8EtFGX6pqDKtvnktr1We3Xgt23E6jK5eL275paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqmoMn3+eda11z7lU5nb8ivNZ3PLLxVl+KWiDL9UlOGXijL8UlGGXyrK8EtFterzR8QF4FXgz8CfMnPx7w3uWN9921n33/a7AvrsKffdr27z2IZcbxX6+PN0cZDP+zPzdx3cj6QVctgvFdU2/Al8NyLORcRWFwVJWo22w/73ZuYrEXEL8L2I+HlmPrP3Bs2bgm8M0shEVydXRMQjwGuZ+bkZt+ntTA53+A2jz8e2zuutT5m50IpZetgfEddHxA1XLwMfBF5c9v4krVabYf8R4NvNu++bgK9n5n93UpWk3nU27F9oYXOG/W2GeWMeIo65tnmGnM+grTGv1z71PuyXNG2GXyrK8EtFGX6pKMMvFWX4paJG1eqrasiv1x5zK6/PadXXuQ1oq0/STIZfKsrwS0UZfqkowy8VZfilogy/VNTaTNE95dNm2/azx9yrH6spv1664pZfKsrwS0UZfqkowy8VZfilogy/VJThl4pamz7/OuvzvPa2hqxtzOtlCtzyS0UZfqkowy8VZfilogy/VJThl4oy/FJRc/v8EXEa+AhwOTPvbK67CfgmcAy4ANyXmb/vr8xds/q2Yz7/uu2542Pu44/1vsduDK/lRbb8XwHuvua6B4HvZ+YdwPebvyVNyNzwZ+YzwJVrrj4BnGkunwHu6bguST1b9jP/kczcAWh+39JdSZJWofdj+yNiC9jqezmSDmfZLf+liDgK0Py+fNANM/NUZm5m5uaSy5LUg2XDfxY42Vw+CTzZTTmSVmXuFN0R8RhwHLgZuAQ8DHwHeBy4HXgZuDczr90puN99tepZjaE9sozKrb4hn7MxT+Hd53pZdIruueHvUtvwr6u2z8GY3/j61OZNdZ3X+aLh9wg/qSjDLxVl+KWiDL9UlOGXijL8UlGT+uruqfb5K+vzORvz8RFT4JZfKsrwS0UZfqkowy8VZfilogy/VJThl4qaVJ9/qtqezz+ktrW3eWxTXm99mrVeNjcX/8Ist/xSUYZfKsrwS0UZfqkowy8VZfilogy/VNRK+/wbGxtsb28f+O9tesZD9oTbnjduP3t/Uz5fv8/ntKvXg1t+qSjDLxVl+KWiDL9UlOGXijL8UlGGXypq7hTdEXEa+AhwOTPvbK57BPgn4LfNzR7KzP+au7CiU3Tbxx9Gn8cBjPk563KK7q8Ad+9z/b9l5l3Nz9zgSxqXueHPzGeAKyuoRdIKtfnM/0BE/CQiTkfEjZ1VJGkllg3/l4B3AncBO8DnD7phRGxFxHZEHHxQv6SVm7vDDyAijgFPXd3ht+i/7XNbd/jtY8w7j6bMHX6zLbXlj4ije/78KPDiMvcjaThzT+mNiMeA48DNEXEReBg4HhF3AQlcAD7ZY42SerDQsL+zha3psL/tOhzzEHLK2jwvU35Oeh32S5o+wy8VZfilogy/VJThl4oy/FJRo5qi2yPhujfldVr1CL1VccsvFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0WNqs8/1d7rmKeSHvM67Xu9TPWxr6put/xSUYZfKsrwS0UZfqkowy8VZfilogy/VNRK+/wbGxtsbx88a1ebfnnbXvuYe/V96vt8/3U9J7/tehvDMQhu+aWiDL9UlOGXijL8UlGGXyrK8EtFGX6pqLlTdEfEbcBXgb8FXgdOZeYXI+Im4JvAMeACcF9m/n7Ofa1ls3zK56WP+fiFKffxh7ToFN2LhP8ocDQzn4+IG4BzwD3Ax4ErmfloRDwI3JiZn5lzX+N9pbVg+Pth+JezaPjnDvszcyczn28uvwqcB24FTgBnmpudYfcNQdJEHOozf0QcA94NPAscycwd2H2DAG7pujhJ/Vn42P6IeCvwBPDpzPzDosOeiNgCtpYrT1Jf5n7mB4iINwNPAU9n5hea634BHM/MnWa/wA8z811z7me8HzBb8DN/P/zMv5zOPvPH7qP8MnD+avAbZ4GTzeWTwJOHLVLScBbZ2/8+4MfAC+y2+gAeYvdz/+PA7cDLwL2ZeWXOfbXazAz5dcdj3kJO1Zi3nlPWWauvS4Zfexn+fnQ27Je0ngy/VJThl4oy/FJRhl8qyvBLRY1qiu55bA2Nj8/JdLnll4oy/FJRhl8qyvBLRRl+qSjDLxVl+KWiJtXnH1KbfvaYTweu2qef8jf1dMUtv1SU4ZeKMvxSUYZfKsrwS0UZfqkowy8VZZ9/BcY880xVY+7jr+oYBLf8UlGGXyrK8EtFGX6pKMMvFWX4paIMv1TU3PBHxG0R8YOIOB8RP42If26ufyQi/i8i/qf5+XD/5fYnM2f+TFVEzPzR+KzqOYsFDig4ChzNzOcj4gbgHHAPcB/wWmZ+buGFRYw2RX65g9ZFZi70Yp17hF9m7gA7zeVXI+I8cGu78iQN7VCf+SPiGPBu4Nnmqgci4icRcToibjzg/2xFxHZEbLeqVFKn5g77/3LDiLcCPwI+m5nfiogjwO+ABP6F3Y8Gn5hzHw77pZ4tOuxfKPwR8WbgKeDpzPzCPv9+DHgqM++ccz+GX+rZouFfZG9/AF8Gzu8NfrMj8KqPAi8etkhJw1lkb//7gB8DLwCvN1c/BNwP3MXusP8C8Mlm5+Cs+5rsln+WtqMCRx3jM+XnpNNhf1cM/3LLHvMLbV1N+TnpbNgvaT0Zfqkowy8VZfilogy/VJThl4pa6Vd3b2xssL3dzyH+bVsvQ7Zu+mwV2oZczro+rr3c8ktFGX6pKMMvFWX4paIMv1SU4ZeKMvxSUas+pfe3wP/uuepmdr8KbIzGWttY6wJrW1aXtf1dZr5tkRuuNPxvWHjEdmZuDlbADGOtbax1gbUta6jaHPZLRRl+qaihw39q4OXPMtbaxloXWNuyBqlt0M/8koYz9JZf0kAGCX9E3B0Rv4iIlyLiwSFqOEhEXIiIF5qZhwedYqyZBu1yRLy457qbIuJ7EfGr5ve+06QNVNsoZm6eMbP0oOtubDNer3zYHxHXAb8EPgBcBJ4D7s/Mn620kANExAVgMzMH7wlHxD8CrwFfvTobUkT8K3AlMx9t3jhvzMzPjKS2RzjkzM091XbQzNIfZ8B11+WM110YYsv/HuClzPx1Zv4R+AZwYoA6Ri8znwGuXHP1CeBMc/kMuy+elTugtlHIzJ3MfL65/CpwdWbpQdfdjLoGMUT4bwV+s+fvi4xryu8EvhsR5yJia+hi9nHk6sxIze9bBq7nWnNnbl6la2aWHs26W2bG664NEf79vh9pTC2H92bmPwAfAj7VDG+1mC8B72R3Grcd4PNDFtPMLP0E8OnM/MOQtey1T12DrLchwn8RuG3P328HXhmgjn1l5ivN78vAt9n9mDIml65Oktr8vjxwPX+RmZcy88+Z+TrwHwy47pqZpZ8AvpaZ32quHnzd7VfXUOttiPA/B9wREe+IiLcAHwPODlDHG0TE9c2OGCLieuCDjG/24bPAyebySeDJAWv5K2OZufmgmaUZeN2NbcbrQQ7yaVoZ/w5cB5zOzM+uvIh9RMTfs7u1h91vNv76kLVFxGPAcXbP+roEPAx8B3gcuB14Gbg3M1e+4+2A2o5zyJmbe6rtoJmln2XAddfljNed1OMRflJNHuEnFWX4paIMv1SU4ZeKMvxSUYZfKsrwS0UZfqmo/wfVHcvvUKWKAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "indices = np.random.choice(binarized_img.size, size=int(binarized_img.size * 0.1), replace=False)\n",
    "noisy_img = np.copy(binarized_img)\n",
    "noisy_img.ravel()[indices] = 1 - noisy_img.ravel()[indices]\n",
    "plt.imshow(noisy_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "markov_random_field = np.array([\n",
    "        [[bn.discrete([0.5, 0.5], name=f\"p(z_({i},{j}))\") for j in range(28)] for i in range(28)], \n",
    "        [[bn.DiscreteVariable(2) for _ in range(28)] for _ in range(28)]])\n",
    "a = 0.9\n",
    "b = 0.9\n",
    "pa = [[a, 1 - a], [1 - a, a]]\n",
    "pb = [[b, 1 - b], [1 - b, b]]\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    bn.discrete(pb, markov_random_field[0, i, j], out=markov_random_field[1, i, j], name=f\"p(x_({i},{j})|z_({i},{j}))\")\n",
    "    if i != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i + 1, j]], name=f\"p(z_({i},{j}), z_({i+1},{j}))\")\n",
    "    if j != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i, j + 1]], name=f\"p(z_({i},{j}), z_({i},{j+1}))\")\n",
    "    markov_random_field[1, i, j].observe(noisy_img[i, j], proprange=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ada61c9f28>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAC4FJREFUeJzt3UHMZWV9x/Hvr6AbZAGhTCcIHWtINy6wEDeaZrrQUDfgAiOrMV2Mi5LUncQNJI2JadS2KxMaiWNSsSSoTEhTJKYtrggDMTJKUWJGHJnM1EwTYWWUfxfvGfI6vO9777z3nnvOO//vJ7m59z1z55z/e977u89zznPufVJVSOrnj6YuQNI0DL/UlOGXmjL8UlOGX2rK8EtNGX6pKcMvNWX4paau3eTGkng5oTSyqsoyz1up5U9yd5JXkrya5MFV1iVps7Lfa/uTXAP8FPgocBZ4Hri/qn6yx/+x5ZdGtomW/0PAq1X186r6LfAt4J4V1idpg1YJ/y3AL7f9fHZY9geSHE9yKsmpFbYlac1WOeG3U9fiHd36qnoEeATs9ktzskrLfxa4ddvP7wVeX60cSZuySvifB25P8r4k7wY+BZxcT1mSxrbvbn9V/S7JA8DTwDXAo1X147VVJmlU+x7q29fGPOaXRreRi3wkHVyGX2rK8EtNGX6pKcMvNWX4paY2+nn+OVs05JksNXoiHRi2/FJThl9qyvBLTRl+qSnDLzVl+KWmHOobOJSnbmz5paYMv9SU4ZeaMvxSU4ZfasrwS00Zfqkpx/k1W37Mely2/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/U1Erj/EnOAG8Avwd+V1V3raMo9bDqDNFeB7CadVzk81dV9es1rEfSBtntl5paNfwFfC/JC0mOr6MgSZuxarf/w1X1epKbgWeS/E9VPbv9CcObgm8M0sxk1ZMub68oeRh4s6q+tMdz1rMxXRXW9drbTdcTflW11C++725/kuuSXH/pMfAx4PR+1ydps1bp9h8CvjO8u14LfLOq/mMtVUka3dq6/UttzG7/gbPJ18e62e3fm0N9UlOGX2rK8EtNGX6pKcMvNWX4pab86u4DYMyPrh7kobxF9vrdug4DbmfLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNOc5/Fbhax+oXjcWv8nv7td+2/FJbhl9qyvBLTRl+qSnDLzVl+KWmDL/UlOP8OrDGvA6gA1t+qSnDLzVl+KWmDL/UlOGXmjL8UlOGX2pqYfiTPJrkQpLT25bdmOSZJD8b7m8Yt8yrW1XteZtSksluGtcyLf/XgbsvW/Yg8P2quh34/vCzpANkYfir6lng4mWL7wFODI9PAPeuuS5JI9vvMf+hqjoHMNzfvL6SJG3C6Nf2JzkOHB97O5KuzH5b/vNJDgMM9xd2e2JVPVJVd1XVXfvclqQR7Df8J4Fjw+NjwJPrKUfSpmSJrzB+DDgK3AScBx4Cvgs8DtwGvAbcV1WXnxTcaV1+xnIHUw/n7eUgD7mtsl8P+O+9VPELw79Ohn9nq/4NDvILdUxjvrbnvM+XDb9X+ElNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS00Zfqkpwy81Zfilpgy/1JThl5pyiu4DYM4fH53SnL8H4SCw5ZeaMvxSU4ZfasrwS00Zfqkpwy81Zfilphzn34Al5kbYUCVXl0X7zesA9mbLLzVl+KWmDL/UlOGXmjL8UlOGX2rK8EtNLQx/kkeTXEhyetuyh5P8KskPh9vHxy1z3qpqz5s0R8u0/F8H7t5h+T9W1R3D7d/XW5aksS0Mf1U9C1zcQC2SNmiVY/4HkvxoOCy4YW0VSdqI/Yb/q8D7gTuAc8CXd3tikuNJTiU5tc9tSRpBljkhleQI8FRVfeBK/m2H516VZ79WPannB3vGMebJ1jn/zapqqeL21fInObztx08Ap3d7rqR5WviR3iSPAUeBm5KcBR4Cjia5AyjgDPCZEWuUNIKluv1r25jd/h3NuQt5kK3ydznIf5NRu/2SDj7DLzVl+KWmDL/UlOGXmjL8UlN+dbdmq+sVeptiyy81Zfilpgy/1JThl5oy/FJThl9qyvBLTTnOr8mM/XFyx/L3ZssvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS005zq9R+Zn8+bLll5oy/FJThl9qyvBLTRl+qSnDLzVl+KWmFo7zJ7kV+AbwJ8BbwCNV9c9JbgT+DTgCnAE+WVX/N16p0xpzvHrRusccz97kFO1XynH8cWWJF95h4HBVvZjkeuAF4F7g08DFqvpikgeBG6rqcwvWNd9X2gJThsTw60pU1VI7bmG3v6rOVdWLw+M3gJeBW4B7gBPD006w9YYg6YC4omP+JEeADwLPAYeq6hxsvUEAN6+7OEnjWfra/iTvAZ4APltVv1m2S5bkOHB8f+VJGsvCY36AJO8CngKerqqvDMteAY5W1bnhvMB/VdWfL1jPfA8wF/CYf/M85t+ftR3zZ+sv8DXg5UvBH5wEjg2PjwFPXmmRkqazzNn+jwA/AF5ia6gP4PNsHfc/DtwGvAbcV1UXF6xrvs3MAnNuIQ8qW/ZxLNvyL9XtXxfDr+0M/zjW1u2XdHUy/FJThl9qyvBLTRl+qSnDLzXlV3drJQ7XHVy2/FJThl9qyvBLTRl+qSnDLzVl+KWmDL/UlOP8S1plPHvOHwd2nL4vW36pKcMvNWX4paYMv9SU4ZeaMvxSU4Zfaspx/g1wLF1zZMsvNWX4paYMv9SU4ZeaMvxSU4ZfasrwS00tDH+SW5P8Z5KXk/w4yd8Nyx9O8qskPxxuHx+/XEnrkkVfNJHkMHC4ql5Mcj3wAnAv8Engzar60tIbS+b7rRbSVaKqlrqqbOEVflV1Djg3PH4jycvALauVJ2lqV3TMn+QI8EHguWHRA0l+lOTRJDfs8n+OJzmV5NRKlUpaq4Xd/refmLwH+G/gC1X17SSHgF8DBfw9W4cGf7NgHXb7pZEt2+1fKvxJ3gU8BTxdVV/Z4d+PAE9V1QcWrMfwSyNbNvzLnO0P8DXg5e3BH04EXvIJ4PSVFilpOsuc7f8I8APgJeCtYfHngfuBO9jq9p8BPjOcHNxrXbNt+ZfYDxuqRFrNWrv962L4pfGtrdsv6epk+KWmDL/UlOGXmjL8UlOGX2pqo+G/8847qapRbqtKsudNutrY8ktNGX6pKcMvNWX4paYMv9SU4ZeaMvxSU5v+SO//Ar/Ytugmtr4KbI7mWttc6wJr26911vanVfXHyzxxo+F/x8aTU1V112QF7GGutc21LrC2/ZqqNrv9UlOGX2pq6vA/MvH29zLX2uZaF1jbfk1S26TH/JKmM3XLL2kik4Q/yd1JXknyapIHp6hhN0nOJHlpmHl40inGhmnQLiQ5vW3ZjUmeSfKz4X7HadImqm0WMzfvMbP0pPtubjNeb7zbn+Qa4KfAR4GzwPPA/VX1k40WsoskZ4C7qmryMeEkfwm8CXzj0mxISf4BuFhVXxzeOG+oqs/NpLaHucKZm0eqbbeZpT/NhPtunTNer8MULf+HgFer6udV9VvgW8A9E9Qxe1X1LHDxssX3ACeGxyfYevFs3C61zUJVnauqF4fHbwCXZpaedN/tUdckpgj/LcAvt/18lnlN+V3A95K8kOT41MXs4NClmZGG+5snrudyC2du3qTLZpaezb7bz4zX6zZF+Hf6Tqw5DTl8uKr+Avhr4G+H7q2W81Xg/WxN43YO+PKUxQwzSz8BfLaqfjNlLdvtUNck+22K8J8Fbt3283uB1yeoY0dV9fpwfwH4DluHKXNy/tIkqcP9hYnreVtVna+q31fVW8C/MOG+G2aWfgL416r69rB48n23U11T7bcpwv88cHuS9yV5N/Ap4OQEdbxDkuuGEzEkuQ74GPObffgkcGx4fAx4csJa/sBcZm7ebWZpJt53c5vxepKLfIahjH8CrgEeraovbLyIHST5M7Zae4BrgW9OWVuSx4CjbH3q6zzwEPBd4HHgNuA14L6q2viJt11qO8oVztw8Um27zSz9HBPuu3XOeL2WerzCT+rJK/ykpgy/1JThl5oy/FJThl9qyvBLTRl+qSnDLzX1/2wNFDzYGLDhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(10000):\n",
    "    i, j = np.random.choice(28, 2)\n",
    "    markov_random_field[1, i, j].send_message(proprange=3)\n",
    "restored_img = np.zeros_like(noisy_img)\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    restored_img[i, j] = np.argmax(markov_random_field[0, i, j].proba)\n",
    "plt.imshow(restored_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
